Load and visualize fits for NEON NO3 models

Intro

This document uses tidybayes to extract and visualize model fit params for the NEON data/ autotrophic uptake project

First created 2/25/2025 by Christa Torrens, with ongoing amendments

Load needed libraries and model fits

# Load libraries
library(tidyverse)
library(tidybayes)
library(bayesplot)
library(rstanarm)
library(here)
library(plotly)
library(sf)

Plotting functions

This keeps all of the core plotting code in 1 place, and allows site-level customization (mostly dataframes and titles)

##### Model fit - Npred + Nobs over time

plot_Nmodelfit_ts <- function(data, start_jday, end_jday, plot_title) {
  plot <- data %>%
    filter(model_jday >= start_jday & model_jday <= end_jday) %>%
    ggplot(aes(x = mod_hours)) +
    geom_point(aes(y = N_conc)) + 
    geom_line(aes(y = N_pred), col = 'red') +
    labs(
      x = "Time (h)",
      y = expression("N" ~ (mmol ~ m^-3))
    ) +
    ggtitle(plot_title) +
    facet_wrap(~yr_jday) +
    theme_bw()
}

plot_Nmodelfit_tsCI <- function(data, start_jday, end_jday, plot_title) {
  plot <- data %>%
    filter(model_jday >= start_jday & model_jday <= end_jday) %>%
    ggplot(aes(x = mod_hours)) +
    geom_ribbon(aes(ymin=N_pred.lower, ymax=N_pred.upper), alpha = 0.7, fill = "slategray") +
    geom_point(aes(y = N_conc)) + 
    geom_line(aes(y = N_pred), col = 'red') +
    labs(
      x = "Time (h)",
      y = expression("N" ~ (mmol ~ m^-3))
    ) +
    ggtitle(plot_title) +
    facet_wrap(~yr_jday) +
    theme_bw()
}

##### Model fit - Npred v Nobs w. 1:1 line

plot_Npred_v_Nobs <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=N_conc, y=N_pred)) +
  geom_point() + 
  labs( x=expression("measured N"~(mmol~m^-3)), # OR (mu*mol~L^-1)
        y=expression("predicted N"~(mmol~m^-3))  # ~ = a space, * = no space
      ) + 
  ggtitle(plot_title) +
  geom_abline(intercept = 0, slope = 1, color = "red") + 
  theme_bw()
}

##### Equifinality checks
# added error bars to all equifinality checks, 4/16/25

## K vs U

plot_KvU <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(y=K, x=U)) +
   geom_errorbar(aes(ymin = K.lower, ymax = K.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
   labs(y=expression("modeled K"~(d^-1)),
       x=expression("modeled U"~(mmol~m^-2~d^-1))
        ) + 
   ggtitle(plot_title) + 
  # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
   theme_bw()

}
  

## K vs N_e

plot_KvNe <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(y=K, x=N_e)) +
   geom_errorbar(aes(ymin = K.lower, ymax = K.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
   labs(y=expression("modeled K"~(d^-1)),
        x=expression("modeled N_e"~(mmol~m^-3)) 
        ) + 
   ggtitle(plot_title) + 
   # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
   theme_bw()
}


## U vs N_e
# Includes CI bars because this relationship tends to look a bit linear - checking whether that holds even w ci values in place

plot_UvNe <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(y=U, x=N_e)) +
   geom_errorbar(aes(ymin = U.lower, ymax = U.upper), 
                width = 0.2, color='grey70') + # CI bars
   geom_point() + 
   labs(y=expression("modeled U"~(mmol~m^-2~d^-1)),
        x=expression("modeled N_e"~(mmol~m^-3)) 
        ) + 
   ggtitle(plot_title) + 
   theme_bw()
}


##### K over time

plot_K_ts <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=d, y=K)) +
   geom_errorbar(aes(ymin = K.lower, ymax = K.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
  # geom_point(y = sumlight.real, color = 'gold') +
   labs( x= "Index", 
         y=expression("modeled K"~(d^-1)) # ~ = a space, * = no space
      ) + 
   ggtitle(plot_title) +
   theme_bw()
}

  
##### U over time

plot_U_ts <- function(data, plot_title) {
  ggplot(data=data, aes(x=d, y=U)) +
   geom_errorbar(aes(ymin = U.lower, ymax = U.upper), 
                 width = 0.2, color='grey70') + # CI bars
   geom_point() + 
  # geom_point(y = sumlight.real_wlou, color = 'gold') +
  # ADD IN HIGH AND LOW CIs
   labs( x= "Index", 
        y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
      ) + 
   ggtitle(plot_title) +
   theme_bw()
}


##### U vs. real sumlight

plot_UvLight <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=sumlightReal, y=U)) +
  geom_point() + 
  labs( x=expression("measured daily light"~(w~m^-2)), 
        y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
      ) + 
  scale_x_log10() + scale_y_log10() + 
  ggtitle(plot_title) +
  theme_bw()
}


##### U + daily (or hourly?) sumlight over time



##### Autotrophic U vs total U

plot_autoUperc <- function(data, plot_title) {
  plot <- ggplot(data=data, aes(x=d, y=U_auto_perc)) +
  geom_point() + 
  xlab("Index") + ylab("Percent autotrophic uptake") + 
  ggtitle(plot_title) + 
  theme_bw()
}

Load model fits and data by site (stored as RDS files), then extract parameters

BIGC

Outputs and visualizations for Upper Big Creek, Fresno, CA

Load model fits and data

########## Load model fit ______________________________________________________

bigc.fit <- readRDS(here("N_uptake_NEON/data/model_fits/bigc.fit.rds"))

# if working from model run
# bigc.fit <- fit.bigc


########## Load model data _____________________________________________________

bigc.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/bigc.data.rds"))


########## Load full daily and hourly datasets _________________________________
# these include datatime info, and the model data does not)

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/bigc_clean.csv")
# bigc.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="US/Pacific"), 
#          model_datetime = local_datetime - hours(4))  


# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/bigc_hourly_clean.csv")
bigc.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Pacific"), 
         model_datetime = local_datetime - hours(4))  

Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a BIGC tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd   2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.05    0.00 0.00   0.05  0.05  0.05  0.05  0.06  2811 1.00
sigma_U  0.16    0.00 0.04   0.09  0.13  0.16  0.18  0.23   129 1.02
b0      -9.15    0.03 0.74 -10.69 -9.62 -9.14 -8.67 -7.77   667 1.01
b1       0.82    0.00 0.08   0.66  0.77  0.82  0.88  1.00   683 1.01

Samples were drawn using NUTS(diag_e) at Tue Apr  1 09:45:41 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Create visualizations

########## N Model fit to data __________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_bigc <- plot_Nmodelfit_ts(data=N_output_bigc.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Upper Big Creek, CA")

# quartz(width=6,height=6)
N_Npred_ts_bigc

### w. credible interval ribbon
N_Npred_tsCI_bigc <- plot_Nmodelfit_tsCI(data=N_output_bigc.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Upper Big Creek, CA")

# quartz(width=6,height=6)
N_Npred_tsCI_bigc

######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_bigc <- plot_Npred_v_Nobs(data=N_output_bigc.df, plot_title="Predicted N vs. measured N, Upper Big Creek, CA")

# quartz(width=6, height=6)
Npred_v_N_bigc

########## Equifinality checks _________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_bigc <- plot_KvU(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_bigc

####### Check equifinality!! K vs N_e  #######  

KvNe_bigc <- plot_KvNe(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_bigc

####### Check equifinality!! U vs N_e  #######  

UvNe_bigc <- plot_UvNe(data=daily_output_bigc.df, plot_title = "Equifinality check, Upper Big Creek, CA: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_bigc

########## Other parameter visualizations _______________________________________

#######  K over time  #######  

K_v_time_bigc <- plot_K_ts(data=daily_output_bigc.df, plot_title = "modeled K over time, Upper Big Creek, CA")

# quartz(width=6.5, height=6)
K_v_time_bigc

##### U over time  #######  

U_time_bigc <- plot_U_ts(data=daily_output_bigc.df, plot_title="Diel nitrate uptake (modeled), Upper Big Creek, CA")

# quartz(width=6, height=6)
U_time_bigc

###### U vs sumlight  #######  

U_vs_light_bigc <- plot_UvLight(data=daily_output_bigc.df, plot_title = "Scatterplot of NO3 uptake and daily light: Upper Big Creek, CA")

# quartz(width=6.5, height=6)
U_vs_light_bigc

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_bigc.df <- daily_output_bigc.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_bigc.df$U_auto_perc) #11.7% on average

Uauto_perc_bigc <- plot_autoUperc(data=daily_output_bigc.df, plot_title = "Daily  autotrophic uptake (as % of total uptake), Upper Big Creek, CA")

# quartz(width=7.5, height=6)
Uauto_perc_bigc 

#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_bigc)

Save updated model output datasets

# path = here("N_uptake_NEON/data/model_output/bigc_output_daily")
# write_csv(daily_output_bigc.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/bigc_output_hourly")
# write_csv(N_output_bigc.df, file=path.h)

CARI

Outputs and visualizations for Caribou Creek, Fairbanks North Star, AK

Load model fits and data

########## Load model fit ______________________________________________________

cari.fit <- readRDS(here("N_uptake_NEON/data/model_fits/cari.fit.rds"))

# if working from model run
# cari.fit <- fit.cari


########## Load model data _____________________________________________________

cari.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cari.data.rds"))


########## Load full hourly datasets __________________________________
# these include datatime info, and the model data does not)

# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/cari_hourly_clean.csv")
cari.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Alaska"), 
         model_datetime = local_datetime - hours(2))# 2-hr offset because it starts getting light at 02:45 in midsummer...

Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a CARI tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=3000; warmup=1500; thin=1; 
post-warmup draws per chain=1500, total post-warmup draws=6000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.10       0 0.00  0.10  0.10  0.10  0.11  0.11 10850    1
sigma_U  0.28       0 0.02  0.23  0.26  0.28  0.29  0.32  1579    1
b0      -2.26       0 0.29 -2.85 -2.45 -2.25 -2.06 -1.69  3517    1
b1       0.19       0 0.03  0.13  0.17  0.19  0.21  0.25  3516    1

Samples were drawn using NUTS(diag_e) at Tue Apr 15 15:00:35 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Create visualizations - Caribou Creek, AK

########## N Model fit to data _________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_cari <- plot_Nmodelfit_ts(data=N_output_cari.df, start_jday = 145, end_jday = 155, plot_title = "NEON N data + model predictions: Caribou Creek, AK")

# quartz(width=6,height=6)
N_Npred_ts_cari

### w. credible interval ribbon
N_Npred_tsCI_cari <- plot_Nmodelfit_tsCI(data=N_output_cari.df, start_jday = 145, end_jday = 155, plot_title = "NEON N data + model predictions: Caribou Creek, AK")


# quartz(width=6,height=6)
N_Npred_tsCI_cari

######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_cari <- plot_Npred_v_Nobs(data=N_output_cari.df, plot_title="Predicted N vs. measured N, Caribou Creek, AK")

# quartz(width=6, height=6)
Npred_v_N_cari

########## Equifinality checks ________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_cari <- plot_KvU(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_cari

####### Check equifinality!! K vs N_e  #######  

KvNe_cari <- plot_KvNe(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_cari

####### Check equifinality!! U vs N_e  #######  

UvNe_cari <- plot_UvNe(data=daily_output_cari.df, plot_title = "Equifinality check, Caribou Creek, AK: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_cari

########## Other parameter visualizations ______________________________________

#######  K over time  #######  

K_v_time_cari <- plot_K_ts(data=daily_output_cari.df, plot_title = "modeled K over time, Caribou Creek, AK")

# quartz(width=6.5, height=6)
K_v_time_cari

##### U over time  #######  

U_time_cari <- plot_U_ts(data=daily_output_cari.df, plot_title="Diel nitrate uptake (modeled), Caribou Creek, AK")

# quartz(width=6, height=6)
U_time_cari

###### U vs sumlight  #######  

U_vs_light_cari <- plot_UvLight(data=daily_output_cari.df, plot_title = "Scatterplot of NO3 uptake and daily light: Caribou Creek, AK")

# quartz(width=6.5, height=6)
U_vs_light_cari

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_cari.df <- daily_output_cari.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_cari.df$U_auto_perc) #4.9% on average

Uauto_perc_cari <- plot_autoUperc(data=daily_output_cari.df, plot_title = "Daily  autotrophic uptake (as % of total uptake), Caribou Creek, AK")

quartz(width=7.5, height=6)
Uauto_perc_cari 

#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_cari)

Save updated model output datasets

# path = here("N_uptake_NEON/data/model_output/cari_output_daily")
# write_csv(daily_output_cari.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/cari_output_hourly")
# write_csv(N_output_cari.df, file=path.h)

CUPE

Outputs and visualizations for Rio Cupeyes, San German Municipio, PR

Load model fits and data

########## Load model fit ______________________________________________________

cupe.fit <- readRDS(here("N_uptake_NEON/data/model_fits/cupe.fit.rds"))

# if working from model run
# cupe.fit <- fit.cupe


########## Load model data _____________________________________________________

cupe.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/cupe.data.rds"))

########## Load full daily and hourly datasets __________________________________
# these include datatime info, and the model data does not)

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/cupe_clean.csv")
# cupe.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="America/Puerto_Rico"), 
#          model_datetime = local_datetime - hours(4))  


# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/cupe_hourly_clean.csv")
cupe.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="America/Puerto_Rico"), 
         model_datetime = local_datetime - hours(4))  

Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a CUPE tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.15    0.00 0.00  0.14  0.15  0.15  0.15  0.15  5813    1
sigma_U  0.24    0.00 0.01  0.22  0.24  0.24  0.25  0.26  2620    1
b0      -5.25    0.01 0.44 -6.08 -5.56 -5.25 -4.96 -4.37  2832    1
b1       0.56    0.00 0.05  0.46  0.53  0.56  0.60  0.66  2841    1

Samples were drawn using NUTS(diag_e) at Tue Apr  1 15:52:19 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Create visualizations - Rio Cupeyes, San German Municipio, PR

########## N Model fit to data __________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_cupe <- plot_Nmodelfit_ts(data=N_output_cupe.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Rio Cupeyes, PR")

# quartz(width=6,height=6)
N_Npred_ts_cupe

### w. credible interval ribbon
N_Npred_tsCI_cupe <- plot_Nmodelfit_tsCI(data=N_output_cupe.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: Rio Cupeyes, PR")

# quartz(width=6,height=6)
N_Npred_tsCI_cupe

######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_cupe <- plot_Npred_v_Nobs(data=N_output_cupe.df, plot_title="Predicted N vs. measured N, Rio Cupeyes, PR")

# quartz(width=6, height=6)
Npred_v_N_cupe

########## Equifinality checks _________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_cupe <- plot_KvU(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_cupe

####### Check equifinality!! K vs N_e  #######  

KvNe_cupe <- plot_KvNe(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_cupe

####### Check equifinality!! U vs N_e  #######  

UvNe_cupe <- plot_UvNe(data=daily_output_cupe.df, plot_title = "Equifinality check, Rio Cupeyes, PR: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_cupe

########## Other parameter visualizations ______________________________________

#######  K over time  #######  

K_v_time_cupe <- plot_K_ts(data=daily_output_cupe.df, plot_title = "modeled K over time, Rio Cupeyes, PR")

# quartz(width=6.5, height=6)
K_v_time_cupe

##### U over time  #######  

U_time_cupe <- plot_U_ts(data=daily_output_cupe.df, plot_title="Diel nitrate uptake (modeled), Rio Cupeyes, PR")

# quartz(width=6, height=6)
U_time_cupe

###### U vs sumlight  #######  

U_vs_light_cupe <- plot_UvLight(data=daily_output_cupe.df, plot_title = "Scatterplot of NO3 uptake and daily light: Rio Cupeyes, PR")

# quartz(width=6.5, height=6)
U_vs_light_cupe

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K * z
# TOTAL UPTAKE = U + (Equilibrium nitrate * K * z)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_cupe.df <- daily_output_cupe.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_cupe.df$U_auto_perc) #4.7% on average

Uauto_perc_cupe <- plot_autoUperc(data=daily_output_cupe.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Rio Cupeyes, PR")

# quartz(width=7.5, height=6)
Uauto_perc_cupe 

#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_cupe)

Save updated model output datasets

# path = here("N_uptake_NEON/data/model_output/cupe_output_daily")
# write_csv(daily_output_cupe.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/cupe_output_hourly")
# write_csv(N_output_cupe.df, file=path.h)

PRIN

Outputs and visualizations for Pringle Creek, Wise County, TX

Load model fits and data

########## Load model fit ______________________________________________________

prin.fit <- readRDS(here("N_uptake_NEON/data/model_fits/prin.fit.rds"))

# if working from model run
# prin.fit <- fit.prin


########## Load model data _____________________________________________________

prin.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/prin.data.rds"))

########## Load full daily and hourly datasets __________________________________
# these include datatime info, and the model data does not)

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/prin_clean.csv")
# prin.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="US/Central"), 
#          model_datetime = local_datetime - hours(4))  


# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/prin_hourly_clean.csv")
prin.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Central"), 
         model_datetime = local_datetime - hours(4))  

Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a PRIN tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.31    0.00 0.00  0.30  0.30  0.31  0.31  0.32  6493    1
sigma_U  0.36    0.00 0.04  0.28  0.33  0.35  0.38  0.44  1104    1
b0      -6.48    0.03 1.26 -9.11 -7.28 -6.48 -5.65 -4.03  1608    1
b1       0.70    0.00 0.15  0.42  0.61  0.70  0.80  1.00  1621    1

Samples were drawn using NUTS(diag_e) at Wed Apr 16 00:03:30 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Create visualizations - Pringle Creek, Wise County, TX

########## N Model fit to data __________________________________________________

######  N model fit over time - predicted N vs measured N  #######

N_Npred_ts_prin <- plot_Nmodelfit_ts(data=N_output_prin.df, start_jday = 90, end_jday = 120, plot_title = "NEON N data + model predictions: Pringle Creek, TX")

# quartz(width=6,height=6)
N_Npred_ts_prin

### w. credible interval ribbon
N_Npred_tsCI_prin <- plot_Nmodelfit_tsCI(data=N_output_prin.df, start_jday = 90, end_jday = 120, plot_title = "NEON N data + model predictions: Pringle Creek, TX")

# quartz(width=6,height=6)
N_Npred_tsCI_prin

######  predicted N vs measured N + 1:1 line  #######

Npred_v_N_prin <- plot_Npred_v_Nobs(data=N_output_prin.df, plot_title="Predicted N vs. measured N, Pringle Creek, TX")

# quartz(width=6, height=6)
Npred_v_N_prin

########## Equifinality checks _________________________________________________

####### Check equifinality!! K vs U  #######  

KvU_prin <- plot_KvU(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled K vs modeled U")
  
# quartz(width=7.5, height=6)
KvU_prin

####### Check equifinality!! K vs N_e  #######  

KvNe_prin <- plot_KvNe(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled K vs modeled equilibrium N")

# quartz(width=7.5, height=6)
KvNe_prin

####### Check equifinality!! U vs N_e  #######  

UvNe_prin <- plot_UvNe(data=daily_output_prin.df, plot_title = "Equifinality check, Pringle Creek, TX: modeled U vs modeled equilibrium N")
 
# quartz(width=7.5, height=6)
UvNe_prin

########## Other parameter visualizations ______________________________________

#######  K over time  #######  

K_v_time_prin <- plot_K_ts(data=daily_output_prin.df, plot_title = "modeled K over time, Pringle Creek, TX")

# quartz(width=6.5, height=6)
K_v_time_prin

##### U over time  #######  

U_time_prin <- plot_U_ts(data=daily_output_prin.df, plot_title="Diel nitrate uptake (modeled), Pringle Creek, TX")

# quartz(width=6, height=6)
U_time_prin

###### U vs sumlight  #######  

U_vs_light_prin <- plot_UvLight(data=daily_output_prin.df, plot_title = "Scatterplot of NO3 uptake and daily light: Pringle Creek, TX")

# quartz(width=6.5, height=6)
U_vs_light_prin

Compare U_total and U_auto

What percentage of total uptake is the autotrophic uptake?

# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K * z
# TOTAL UPTAKE = U + (Equilibrium nitrate * K * z)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_prin.df <- daily_output_prin.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)

U_auto_avg <- mean(daily_output_prin.df$U_auto_perc) #12.74% on average

Uauto_perc_prin <- plot_autoUperc(data=daily_output_prin.df, plot_title = "Daily autotrophic uptake (as % of total uptake), Pringle Creek, TX")

# quartz(width=7.5, height=6)
Uauto_perc_prin 

#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_prin)

Save updated model output datasets

# path = here("N_uptake_NEON/data/model_output/prin_output_daily")
# write_csv(daily_output_prin.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/prin_output_hourly")
# write_csv(N_output_prin.df, file=path.h)

WLOU

Ouptuts and visualizations for West St Louis Creek, Grand, CO

Load model fits and data

########## Load model fit ______________________________________________________

wlou.fit <- readRDS(here("N_uptake_NEON/data/model_fits/wlou.fit.rds"))

# if working from model run
# wlou.fit <- fit.wlou


########## Load model data _____________________________________________________

# data from model

wlou.modeldata <- readRDS(here("N_uptake_NEON/data/model_datalist/wlou.data.rds"))

###### Load dataset (includes datatime info, which model data does not) _______

# # 15-min dataset
# path <- here("N_uptake_NEON/data/neon_data_clean/wlou_clean.csv")
# wlou.df <- read_csv(path) %>%
#   mutate(local_datetime = with_tz(local_datetime, tzone="US/Mountain"), 
#          model_datetime = local_datetime - hours(4))  
# 

# hourly dataset
path_h <- here("N_uptake_NEON/data/neon_data_clean/wlou_hourly_clean.csv")
wlou.df.h <- read_csv(path_h) %>%
  mutate(local_datetime = with_tz(local_datetime, tzone="US/Mountain"), 
         model_datetime = local_datetime - hours(4))  

Extract model parameters and real-data values

First, check the fit for sigmas and betas (“sigma”, “sigma_U”, “b0”, “b1”), since these are notoriously tough to fit. Then put the spread and mean/ median/ quartile info for K, U, and N_e into a WLOU tibble (to be saved in one large csv at the end of the document). Finally, assemble tibbles for the visualizations, which include both modeled parameters and real data.

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

         mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma    0.05    0.00 0.00  0.05  0.05  0.05  0.05  0.05  5201    1
sigma_U  0.41    0.00 0.03  0.35  0.39  0.41  0.43  0.47  2253    1
b0      -2.64    0.01 0.57 -3.74 -3.03 -2.64 -2.25 -1.52  3282    1
b1       0.11    0.00 0.07 -0.02  0.07  0.11  0.16  0.24  3283    1

Samples were drawn using NUTS(diag_e) at Mon Mar 10 15:19:42 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Create visualizations - West St. Louis Creek, CO

########## N Model fit to data ________________________________________________

######  N model fit over time  #######

# N_and_Npred_wlou <- N_output_wlou.df %>%
#   #group_by(year(model_datetime_wlou)) %>%
#   filter(model_jday >= 35 & model_jday <= 60) %>%  # to see the little boxes better...
#   ggplot(aes(x=mod_hours)) +
#   geom_point(aes(y=N_conc)) + 
#   geom_line(aes(y=N_conc_pred), col='red')+
#   labs(
#     x="Time (h)", y=expression("N"~(mmol~m^-3))
#   ) +
#   ggtitle("NEON: West St. Louis Creek N data + model predictions") +
#   facet_wrap(~yr_jday) +
#   theme_bw()

N_Npred_ts_wlou <- plot_Nmodelfit_ts(data=N_output_wlou.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: West St. Louis Creek, CO ")

#quartz(width=6,height=6)
N_Npred_ts_wlou

### w. credible interval ribbon
N_Npred_tsCI_wlou <- plot_Nmodelfit_tsCI(data=N_output_wlou.df, start_jday = 35, end_jday = 60, plot_title = "NEON N data + model predictions: West St. Louis Creek, CO")

# quartz(width=6,height=6)
N_Npred_tsCI_wlou

######  N-pred vs N  #######

# Npred_v_N_wlou <- ggplot(data = N_output_wlou.df, aes(x=N_conc, y=N_conc_pred)) +
#   geom_point() + 
#   labs( x=expression("measured N"~(mu*mol~L^-1)), # fcn changed both to mmol~m^-3
#         y=expression("modeled N"~(mu*mol~L^-1))  # ~ = a space, * = no space
#       ) + 
#   ggtitle("Measured N vs predicted N, West St. Louis Creek, CO") +
#   geom_abline(intercept = 0, slope = 1, color = "red") + 
#   theme_bw()

Npred_v_N_wlou <- plot_Npred_v_Nobs(data=N_output_wlou.df, plot_title="Predicted N vs. measured N, West St. Louis Creek, CO")

#quartz(width=6, height=6)
Npred_v_N_wlou

########## Equifinality checks ________________________________________________

####### Check equifinality!! K vs U  #######  

# KvU_wlou <- ggplot(data=daily_output_wlou.df, aes(y=K_mod_avg_wlou, x=U_mod_avg_wlou)) +
#   geom_point() + 
#   labs(y=expression("modeled K"~(d^-1)),
#        x=expression("modeled U"~(mmol~m^-2~d^-1))
#         ) + 
#   ggtitle("Equifinality check, West St. Louis Creek, CO: modeled K vs modeled U") + 
#   # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
#   theme_bw()


KvU_wlou <- plot_KvU(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled K vs modeled U")

  
#quartz(width=7.5, height=6)
KvU_wlou

####### Check equifinality!! K vs N_e  #######  

# KvNe_wlou <- ggplot(data=daily_output_wlou.df, aes(y=K, x=N_e)) +
#    geom_point() + 
#    labs(y=expression("modeled K"~(d^-1)),
#         x=expression("modeled N_e"~(mmol~m^-3)) 
#         ) + 
#    ggtitle("Equifinality check, West St. Louis Creek, CO: modeled K vs modeled equilibrium N") + 
#    # geom_abline(intercept = 0, slope = 1, color = "red") + # doesn't show up - off the charts!
#   theme_bw()


KvNe_wlou <- plot_KvNe(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled K vs modeled equilibrium N")

#quartz(width=7.5, height=6)
KvNe_wlou

####### Check equifinality!! U vs N_e  #######  

# UvNe_wlou <- ggplot(data=daily_output_wlou.df, aes(y=U, x=N_e)) +
#    geom_point() + 
#    labs(y=expression("modeled U"~(mmol~m^-2~d^-1)),
#         x=expression("modeled N_e"~(mmol~m^-3)) 
#         ) + 
#    ggtitle("Equifinality check, West St. Louis Creek, CO: modeled U vs modeled equilibrium N") + 
#   theme_bw()

UvNe_wlou <- plot_UvNe(data=daily_output_wlou.df, plot_title = "Equifinality check, West St. Louis Creek, CO: modeled U vs modeled equilibrium N")
 
#quartz(width=7.5, height=6)
UvNe_wlou

########## Other parameter visualizations ______________________________________

######## K over time  #######  

# K_v_time_wlou <- ggplot(data = daily_output_wlou.df, aes(x=d, y=K)) +
#   geom_point() + 
#   # geom_point(y = sumlight.real, color = 'gold') +
#   # ADD IN HIGH AND LOW CIs
#   # xlab("Julian day") + ylab("modeled K (day -1)") + # daily change in N concentration
#   labs( x= "Index", 
#         y=expression("modeled K"~(d^-1)) # ~ = a space, * = no space
#       ) + 
#   #ylim(0,20) +
#   ggtitle("modeled K over time, West St. Louis Creek, CO") +
#   theme_bw()


K_v_time_wlou <- plot_K_ts(data=daily_output_wlou.df, plot_title = "modeled K over time, West St. Louis Creek, CO")

#quartz(width=6.5, height=6)
K_v_time_wlou

##### to clip: no function yet
# K_time_clip_wlou <- ggplot(data = daily_output_wlou.df, aes(x=mod_day_wlou, y=K_mod_avg_wlou)) +
#   geom_point() + 
#   # geom_point(y = sumlight.real, color = 'gold') +
#   # ADD IN HIGH AND LOW CIs
#   xlab("Julian day") + ylab("modeled K (day -1)") + # ??? UNITS?
#   ylim(0,10) +
#   ggtitle("modeled K over time, West St. Louis Creek, CO - ylim = 0,10") +
#   theme_bw()
# 
# quartz()
# K_time_clip_wlou


######## U over time  #######  

# U_time_wlou <- ggplot(data = daily_output_wlou.df, aes(x=d, y=U)) +
#   geom_point() + 
#   # geom_point(y = sumlight.real_wlou, color = 'gold') +
#   # ADD IN HIGH AND LOW CIs
#   #xlab("Julian day") + ylab("modeled U (mmol/m2/day)") + 
#   labs( x= "Index", 
#         y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
#       ) + 
#   #ylim(0,1) +
#   #ggtitle("modeled U over time, Big Creek 2019 pooled model w real light") +
#   ggtitle("Diel nitrate uptake (modeled), West St. Louis Creek, CO") +
#   theme_bw()

U_time_wlou <- plot_U_ts(data=daily_output_wlou.df, plot_title="Diel nitrate uptake (modeled), West St. Louis Creek, CO")

#quartz(width=6, height=6)
U_time_wlou

# U_time_clip_wlou <- ggplot(data = daily_output_wlou.df, aes(x=mod_day_wlou, y=U_mod_avg_wlou)) +
#    geom_point() +
#    # geom_point(y = sumlight.real, color = 'gold') +
#    # ADD IN HIGH AND LOW CIs
#    xlab("Julian day") + ylab("modeled U (mmol/m2/day)") +
#    ylim(0,1) +
#    ggtitle("modeled U over time, West St. Louis Creek, CO - ylim = 0-1") +
#    theme_bw()
# 
#  quartz()
#  U_time_clip_wlou



######## U vs sumlight  #######  

# U_vs_light_wlou <- ggplot(data = daily_output_wlou.df, aes(x=sumlightReal, y=U)) +
#   geom_point() + 
#   #xlab("true light (satellite)") + ylab("modeled U (mmol/m2/day)") +
#   # xlab("light (satellite)") + ylab("modeled U (mmol/m2/day)") +
#   labs( x= "light (satellite)", 
#         y=expression("modeled U"~(mmol~m^-2~d^-1)) # ~ = a space, * = no space
#       ) + 
#   scale_x_log10() + scale_y_log10() + 
#   ggtitle("Scatterplot of NO3 uptake and daily light: West St. Louis Creek, CO") +
#   #ylim = c(-0.2, 1) +
#   theme_bw()

U_vs_light_wlou <- plot_UvLight(data=daily_output_wlou.df, plot_title = "Scatterplot of NO3 uptake and daily light: West St. Louis Creek, CO")

# quartz(width=6.5, height=6)
U_vs_light_wlou

Compare U_total and U_auto

# NON-AUTOTROPHIC UPTAKE = Equilibrium nitrate * K
# TOTAL UPTAKE = U + (Equilibrium nitrate * K)  
# Since U_auto is very small, there won't be much difference between 'non-autotrophic uptake' and 'total uptake'

daily_output_wlou.df <- daily_output_wlou.df %>%
  mutate(U_het = N_e*K*z,      
         U_tot = U_het + U, 
         U_auto_perc = 100*U/U_tot)


U_auto_avg <- mean(daily_output_wlou.df$U_auto_perc) #10.3% on average

# Uauto_plot_wlou <- ggplot(data=daily_output_wlou.df, aes(x=d, y=U_auto_perc)) +
#   geom_point() + 
#   xlab("Index") + ylab("Percent autotrophic uptake") + 
#   ggtitle("Percent of autotrophic uptake each day, West St. Louis Creek, CO") + 
#   theme_bw()


Uauto_perc_wlou <- plot_autoUperc(data=daily_output_wlou.df, plot_title = "Daily  autotrophic uptake (as % of total uptake), West St. Louis Creek, CO")

# quartz(width=7.5, height=6)
Uauto_perc_wlou 

#### Interactive w. ggplotly, but only gives the x,y info... maybe if I added datetime-associated labels?
# ggplotly(Uauto_plot_wlou)

Save updated model output datasets

# path = here("N_uptake_NEON/data/model_output/wlou_output_daily")
# write_csv(daily_output_wlou.df, file=path)
# 
# path.h = here("N_uptake_NEON/data/model_output/wlou_output_hourly")
# write_csv(N_output_wlou.df, file=path.h)

Assemble and save mean, median, 95% QI tibble for all sites

Combine the site tibbles into one large tibble, then save it

########## Reload all datasets as needed  _______________________________________


########## Bind then save daily parameters (K, U, N_e)  _________________________

# daily_summary_all.df <- bind_rows()
# path <- here()
# saveRDS(daily_summary_all.df, path)

########## Bind then save hourly parameters (N_pred)  ___________________________

# hourly_summary_all.df <- bind_rows()
# path2 <- here()
# saveRDS(hourly_summary_all.df, path2)

Create whole-dataset visualizations

########## Violin plots _______________________________________________

##### K


##### U


##### N_e